/*
   FieldVis - A Tool for Visualizing Astrophysical Magnetohydrodynamic Flow
   
   Copyright (C) 2006 - 2008     FieldVis Development Team
  
   This program is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 2 of the License, or
   (at your option) any later version.
  
   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.
  
   You should have received a copy of the GNU General Public License along
   with this program; if not, write to the Free Software Foundation, Inc.,
   51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

   File     : Streamtube.cpp
   Author   : Blayne Field
   
   Purpose  : Implementation file for the Streamtube class.  This class holds
              all of the relevant information for creating a streamtube.
 */

#include "Streamtube.h"


Streamtube::Streamtube(VectorField *vfield, Point seed, int maxPoints, bool BothDirections,
                       float minThreshold, int ID, ColorMap *pColorMap, int NumSides, Texture *pTex,
                       bool proportionalTexture, ScalarField *helicity, float Radius)
   : vectorField(vfield), MaxIntegrationPoints(maxPoints), bothDirections(BothDirections),
     MinimumThreshold(minThreshold), pointArraySize(0), pointArrayAllocated(0),
     indexArraySize(0), indexArrayAllocated(0), sides(NumSides), id(ID), radius(Radius),
     colorMap(pColorMap), texture(pTex), helicityField(helicity), normalArray(0),
     vertexArray(0), colorArray(0), textureArray(0), indices(0)
{
   PutTubeAt(seed);
}

// destructor
Streamtube::~Streamtube()
{
   if (indices)
      delete [] indices;

   if (normalArray)
      delete [] normalArray;
}

void Streamtube::PutTubeAt(Point seedPoint)
{
   sline = ComputeTubePath(seedPoint);

   points = sline.GetPointsCalculated();

   if (points > 1)
   {
      AllocatePointArray();

      ComputeNormalArray();  
      ComputeVertexArray();
      ComputeColorArray();
      ComputeTextureArray(vectorField);

      AllocateIndexArray();

      ComputeIndexArray();
   }
}

// return a path based on a given seed point
StreamlinePath Streamtube::ComputeTubePath(Point seedPoint)
{
   Integrator rk(vectorField);
   StreamlinePath newPath;
   
   if (bothDirections)
   {
      newPath = rk.CalculateBidirectionalPath(seedPoint, MaxIntegrationPoints, MinimumThreshold);
   }
   else
   {
      newPath = rk.CalculatePath(seedPoint, MaxIntegrationPoints, FORWARD, MinimumThreshold);
   }

   newPath.ScaleByFactor(1.0 / vectorField->GetScaleFactor());
   newPath.SmoothPath(vectorField);

   return (newPath);
}

void Streamtube::AllocatePointArray()
{
// Combine the different arrays into one...
// localize the call to new here for the GLfloat arrays....
//   normalArray = new GLfloat[points*(sides+1)*3];
//   vertexArray = new GLfloat[points*(sides+1)*3];
//   colorArray = new GLfloat[points*(sides+1)*4];
//   textureArray = new GLfloat[points*(sides+1)*2];
//  +---------------|---------------|--------------------|----------+
//    normalArray     vertexArray      colorArray          textureArray
//  +---------------|---------------|--------------------|----------+

   pointArraySize = (points * (sides + 1)) * 12;
      
   if (pointArraySize > pointArrayAllocated)
   {
      pointArrayAllocated = (pointArraySize * 5) / 4;    // 25% bigger
         
      if (normalArray)
         delete [] normalArray;
          
      normalArray = new GLfloat[ pointArrayAllocated ];      
   }
      
// set pointers to the start of the internal arrays
   vertexArray = normalArray + (points * (sides + 1)) * 3;
   colorArray = vertexArray + (points * (sides + 1)) * 3;
   textureArray = colorArray + (points * (sides + 1)) * 4;
}

void Streamtube::AllocateIndexArray()
{
   indexArraySize = (points - 1) * (sides + 1) * 2;
   
   if (indexArraySize > indexArrayAllocated)
   {
      indexArrayAllocated = (indexArraySize * 5) / 4;    // 25% bigger

      if (indices)
         delete [] indices;

      indices = new GLuint[ indexArrayAllocated ];      
   }
}

void Streamtube::ComputeIndexArray()
{
   GLuint *pIndex = indices;

   for (int i = 0; i < (points-1)*(sides+1); i++)
   {
      *pIndex++ = i;
      *pIndex++ = i + sides + 1;
   }
}

void Streamtube::ComputeColorArray()
{
   Point *slCoords = sline.GetCoordinates();
   Vector3D tempColor;
   GLfloat *pColor = colorArray;

   for (int point = 0; point < points; point++)
   {
      tempColor = colorMap->GetRGBValue(slCoords[point]);
      for (int side = 0; side <= sides; side++)
      {
         *pColor++ = tempColor.GetX();
         *pColor++ = tempColor.GetY();
         *pColor++ = tempColor.GetZ();      
         *pColor++ = 1.0;
      }
   }
}

void Streamtube::DrawStreamtube(GLenum drawMode)
{
   if (points > 1)
   {
      glPushAttrib(GL_ENABLE_BIT);
      glEnable(GL_TEXTURE_2D);
      glEnableClientState(GL_VERTEX_ARRAY);
      glVertexPointer(3, GL_FLOAT, 0, vertexArray);
      glEnableClientState(GL_COLOR_ARRAY);
      glColorPointer(4, GL_FLOAT, 0, colorArray);
      glEnableClientState(GL_NORMAL_ARRAY);
      glNormalPointer(GL_FLOAT, 0, normalArray);
      glEnableClientState(GL_TEXTURE_COORD_ARRAY);
      glTexCoordPointer(2, GL_FLOAT, 0, textureArray);

      if (UseTexture())
      {
         glBindTexture(GL_TEXTURE_2D, texture->GetID());
      }

      glPolygonMode(GL_FRONT_AND_BACK, drawMode);
      glPushName(id);
      for (volatile int i = 0; i < points - 1; i++)
      {
         glPushName(i);
         glDrawElements(GL_QUAD_STRIP, (sides+1)*2, GL_UNSIGNED_INT, indices + (i*(sides+1)*2));
         glPopName();     
      }
      glPopName();
      glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);

      glDisableClientState(GL_VERTEX_ARRAY);
      glDisableClientState(GL_COLOR_ARRAY);
      glDisableClientState(GL_NORMAL_ARRAY);
      glDisableClientState(GL_TEXTURE_COORD_ARRAY);
      glDisable(GL_TEXTURE_2D);

      // draw the cap
      if (sline.GetExitsGrid())
      {
         Vector3D lastDir = sline.GetLastDirection();
         lastDir.Normalize();
         lastDir = (lastDir * DEFAULT_CAP_LENGTH);
         Point tip = sline.GetLastCoordinate() + lastDir;
         lastDir.Normalize();

         // push the id and last side onto the name stack - this prevents a program
         // crash when the user clicks on a cap.  now the last point before the cap is
         // highlighted instead.
         glPushName(id);
         glPushName(points - 1);
         for (int side = 0; side <= sides; side++)
         {
            //points*(sides+1)*3
            Point vertex1 = sline.GetLastCoordinate();
            int index = ArrayIndex(points-1, side);
            vertex1.MoveX(normalArray[index]*2*radius);
            vertex1.MoveY(normalArray[index+1]*2*radius);
            vertex1.MoveZ(normalArray[index+2]*2*radius);

            Point vertex2 = sline.GetLastCoordinate();
            index = ArrayIndex(points-1, (side+1)%sides);
            vertex2.MoveX(normalArray[index]*2*radius);
            vertex2.MoveY(normalArray[index+1]*2*radius);
            vertex2.MoveZ(normalArray[index+2]*2*radius);

            Vector3D tipNormal = 
            Point::VectorFromTo(vertex1, tip).Cross(Point::VectorFromTo(vertex1, vertex2));
            tipNormal.Normalize();

            glBegin(GL_POLYGON);
            glColor3f(1,0,0);
            glVertex3d(tip.GetX(), tip.GetY(), tip.GetZ());
            glNormal3d(tipNormal.GetX(), tipNormal.GetY(), tipNormal.GetZ());
            glVertex3d(vertex1.GetX(), vertex1.GetY(), vertex1.GetZ());
            glNormal3d(tipNormal.GetX(), tipNormal.GetY(), tipNormal.GetZ());
            glVertex3d(vertex2.GetX(), vertex2.GetY(), vertex2.GetZ());
            glNormal3d(tipNormal.GetX(), tipNormal.GetY(), tipNormal.GetZ());
            glEnd();

         }
         glPopName();         
         glPopName();         
      }
      glPopAttrib();
   }
}

void Streamtube::ComputeNormalArray()
{
   // void Streamtube::ComputeVectorsForRotating()

   Vector3D *startVectors = new Vector3D[2 * points];
   Vector3D *planeNormals = startVectors + points;

   Vector3D *directions = sline.GetDirections();

   // the start vectors are where the rotation around the axis starts
   // the first one just needs to be perpendicular to the initial direction
   startVectors[0].SetX(directions[0].GetY());
   startVectors[0].SetY(-(directions[0].GetX()));
   startVectors[0].SetZ(0);
   startVectors[0].Normalize();

   // for the remaining start vectors, find one in the plane formed by the
   // axis of the tube and the previous start vector
   for (int i = 1; i < points; i++)
   {
      // get the normal vector for the previous plane 
      planeNormals[i-1] = directions[i-1].Cross(startVectors[i-1]);
      planeNormals[i-1].Normalize();

      // get a vector perpendicular to the normal of the previous plane
      // and the current direction
      startVectors[i] = planeNormals[i-1].Cross(directions[i]);
      startVectors[i].Normalize();
   } 

   // the last segment of the tube needs to be set
   planeNormals[points - 1] = planeNormals[points - 2];
   planeNormals[points - 1].Normalize();


   // void Streamtube::ComputeNormalArray()
   // calculate the necessary trig functions only once
   // add extra element at the end because the fill loop below goes to <= sides
   // and we don't need the extra modulus operation......
   float *cosArray = new float[2*(sides+1)];
   float *sinArray = cosArray + (sides + 1);
   for (int i = 0; i < sides; i++)
   {
      cosArray[i] = cos(i * 2 * M_PI / sides);
      sinArray[i] = sin(i * 2 * M_PI / sides);
   }
   cosArray[sides] = cosArray[0];
   sinArray[sides] = sinArray[0];

   GLfloat *pNorm = normalArray;
   Vector3D tempNormal;

   for (int point = 0; point < points; point++)
   {
      // these two vectors are perpendicular to each other
      Vector3D v = startVectors[point];
      Vector3D u = planeNormals[point];     

      for (int side = 0; side <= sides; side++)
      {
         tempNormal = (u * cosArray[side]) + (v * sinArray[side]);
         tempNormal.Normalize();

         *pNorm++ = tempNormal.GetX();
         *pNorm++ = tempNormal.GetY();
         *pNorm++ = tempNormal.GetZ();     
      }
   }

   delete [] cosArray;
   delete [] startVectors;
}

void Streamtube::ComputeVertexArray()
{
   Point *slCoords = sline.GetCoordinates();

   GLfloat *pVert = vertexArray;
   GLfloat *pNorm = normalArray;

   for (int point = 0; point < points; point++)
   {
      Point currentPoint = slCoords[point];

      for (int side = 0; side <= sides; side++)
      {
         *pVert++ = currentPoint.GetX() + *pNorm++ * radius;
         *pVert++ = currentPoint.GetY() + *pNorm++ * radius;
         *pVert++ = currentPoint.GetZ() + *pNorm++ * radius;
      }
   }
}

void Streamtube::ComputeTextureArray(VectorField *stretch)
{
   //void Streamtube::ComputeHelicity()

   float *textureStart = new float[2 * points];
   float *textureYStart = textureStart + points;

   Point *slCoords = sline.GetCoordinates();

   if ( ! UseHelicity())
   {
      for (int i = 0; i < points; i++)
         textureStart[i] = 0.0;
   }
   else
   {
      float dphi;
      textureStart[0] = 0.0;
      for (int i = 1; i < points; i++)
      {
         float factor = helicityField->GetScaleFactor();
         float dist = factor * Point::DistanceBetween(slCoords[i], slCoords[i-1]);
         float val1 = helicityField->GetInterpolatedValue(slCoords[i] * factor);
         float val2 = helicityField->GetInterpolatedValue(slCoords[i] * factor);

         dphi = 0.5 * dist * (val1 + val2);
         textureStart[i] = textureStart[i-1] + (dphi / (2 * M_PI)) * sides;
      }
   }

   // void Streamtube::ComputeTextureStretch(VectorField *stretch)

   for (int i = 0; i < points; i++)
   {
      textureYStart[i] = 0.0;
   }
   float stretchFactor;
   bool shouldStretch = true;

   // no point in calculating the amount to stretch if there's no texture
   if (UseTexture())
   {
      for (int i = 1; i < points; i++)
      {
         float mag = stretch->GetMagnitude(slCoords[i-1] * stretch->GetScaleFactor());
         float dist = stretch->GetScaleFactor() * 
                      Point::DistanceBetween(slCoords[i], slCoords[i-1]);

         float minDist = 1 / texture->GetMinMaxFactor();
         if (shouldStretch)
         {
            if (mag > texture->GetMaxValue())
               stretchFactor = 1;
            else if (mag < texture->GetMinValue())
               stretchFactor = minDist;
            else
               stretchFactor = minDist + (1 - minDist) * (mag - texture->GetMinValue()) /
                               (texture->GetMaxValue() - texture->GetMinValue());
         }
         else
            stretchFactor = 1;

         textureYStart[i] = textureYStart[i-1] + dist /
                            (texture->GetTextureLength() * stretchFactor * stretch->GetScaleFactor());
      }  
   }

   // void Streamtube::ComputeTextureArray()

   GLfloat *pTexture = textureArray;

   for (int point = 0; point < points; point++)
   {
      for (int side = 0; side <= sides; side++)
      {
         *pTexture++ = side + textureStart[point];
         *pTexture++ = textureYStart[point];
      }
   }

   delete [] textureStart;
}

